Function-on-Scalar regressions for the white mold data

Author

Denis Shah

Published

March 4, 2025

Packages

library(fda)
library(refund)
library(fdatest)
library(tidyverse)
library(tictoc)
library(cowplot)
library(kableExtra)

Load and process the data

wm_load <- readr::read_csv(here::here("Data", "data_model_plus_weather_filtered.csv"), show_col_types = FALSE)

# Simplify things. Keep only the weather-related variables and others needed for calculations
wm_data <-
  wm_load %>%
  dplyr::select(subject, date, planting.date, sampling.date, wm, d2m, t2m_mean, t2m_max, t2m_min, st, sp, sm, rh) %>%
  # Do the filtering steps before doing any calculations or feature engineering:
  # wm_load has identical data for each of the sampling dates, which is why there is a filtering step on sampling.date.
  dplyr::group_by(subject) %>% 
  dplyr::filter(sampling.date == max(sampling.date)) %>%
  dplyr::ungroup() %>% 
  # Add the response variable (wm present or absent; binary):
  dplyr::group_by(subject) %>% 
  dplyr::mutate(wm = (mean(wm, na.rm = T) > 0)*1) %>% 
  # wm as a factor:
  dplyr::mutate(wm = factor(wm, levels = c(0, 1))) %>%
  dplyr::ungroup() %>% 
  dplyr::filter(!is.na(wm)) %>% 
  # Calculate dap (as a numeric):
  dplyr::mutate(dap = as.numeric(date - planting.date)) %>%
  # Convert temperatures from Kelvin to degrees Celsius:
  dplyr::mutate(across(d2m:st, ~ .x - 273.15)) %>%
  # dewpoint depression:
  dplyr::mutate(dpd = t2m_mean - d2m) %>%
  # surface pressure in kPa:
  dplyr::mutate(sp = sp/1000) %>%
  # Ratio of soil temperature to soil moisture:
  dplyr::mutate(stsm = st/sm) %>%
  # Difference in max and min temperatures:
  dplyr::mutate(tdiff = t2m_max - t2m_min) %>%
  # Calculate saturation vapor pressure (es):
  dplyr::mutate(es = 0.61078 * exp((17.27 * t2m_mean) / (t2m_mean + 237.3))) %>% 
  # Calculate actual vapor pressure (ea):
  dplyr::mutate(ea = (rh / 100) * es) %>% 
  # Calculate VPD (kPa):
  dplyr::mutate(vpd = es - ea) %>% 
  # estimate GDD (base 0). NB: base 0 is reasonable for snap bean; see the Jenni et al. (2000) paper. I think we want GDD to start accumulating the day after planting date onwards...
  # I use the day after planting, because we don't know exactly what time of the day the field was planted.
  dplyr::mutate(gddi = ifelse(dap <= 0, 0, (t2m_max + t2m_min)*0.5 - 0)) %>%
  dplyr::group_by(subject) %>% 
  # We don't want the gddi column after creating gdd:
  dplyr::mutate(gdd = cumsum(gddi), .keep = "unused") %>%
  dplyr::ungroup() %>%
  # Keep only the columns you need:
  dplyr::select(subject, dap, wm:gdd, -es, -ea) %>%
  # Filter the dataset to dap <= 50:
  dplyr::filter(dap <= 50) %>%
  # wm as a contrast: -1 = absent, 1 = present (needed for function-on-scalar regression):
  dplyr::mutate(wm = as.numeric(wm)-1) %>%
  dplyr::mutate(wm = 2*wm-1)

data.prep <- function(x) {
  # Prepare the data for function-on-scalar regression
  # N x J matrix of the functional predictor 
  # Args:
  #  x = unquoted variable name
  # Returns:
  #  a list with the following named elements:
  #   y = N x J matrix of the functional response
  #   x = a vector of the wm status coded as -1, 1
  #   yind =  a vector of the evaluation days, which is of length 81 from -30 to 50
  
  .x <- enquo(x)
  
  # The environmental data:
  ez <-
    wm_data %>%
    dplyr::select(subject, dap, !!.x) %>%
    tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = !!.x) %>%
    dplyr::select(-subject) %>%
    as.matrix()
  
  colnames(ez) <- NULL
  
  # the vector of wm status:
  x.wm <- 
    wm_data %>%
    dplyr::group_by(subject) %>%
    dplyr::summarise(wm = mean(wm)) %>%
    dplyr::pull(wm)
  
  # the vector of the evaluation days:
  days <- seq(-30, 50)
  
  # the final data frame:
  df <- list(y = ez, x = x.wm, yind = days)
  
  return(df)
  
}

# Example of use:
# dat <- data.prep(d2m)

fANOVA <- function(dat) {
  # Performs a functional ANOVA and plots the estimated coefs for the overall mean and 
  # the difference between wm(0) and wm(1).
  #
  # Args:
  #   dat: prepared data from calling the data.prep function
  #
  # Returns:
  #   a ggplot graphic of the estimated beta(t) coefs.
  
  
  # k = 30 gives wiggliness, no oversmoothing
  m2 <- pffr(y~x, yind = dat$yind, data = dat,
             bs.yindex = list(bs = "ps", k = 30, m = c(2, 1)), 
             bs.int = list(bs = "ps", k = 30, m = c(2, 1)))
  
  
  # The smooth coefficients for the overall mean and beta(t):
  z <- coef(m2)$smterms
  
  # 1. For the overall mean:
  a1 <- z[[1]]$coef[, "yindex.vec"]
  a2 <- z[[1]]$coef[, "value"]
  a3 <- z[[1]]$coef[, "se"]
  
  # 2. For beta(t):
  b1 <- z[[2]]$coef[, "yindex.vec"]
  b2 <- z[[2]]$coef[, "value"]
  b3 <- z[[2]]$coef[, "se"]

  
  # Extract parts of the output to produce a more attractive plot in ggplot.
  # The coefs for the overall mean, adding an approx. 95% CI:
  X <- data.frame(x = a1, coef.mean = a2, lower.mean = a2 - 1.96*a3, upper.mean = a2 + 1.96*a3, coef = "Overall mean")
  
  # The coefs for the time trend of the mean difference between epidemics and non-epidemics:
  Y <- data.frame(x = b1, coef.mean = b2, lower.mean = b2 - 1.96*b3, upper.mean = b2 + 1.96*b3, coef = "Difference")
  
  Z <- rbind(X, Y)
  
  # To get different lines in each facet, you need another data.frame:
  hline.data <- data.frame(z = c(0), coef = c("Difference"))
  
  breaks <- seq(-30, 50, 10)
  labels <- seq(-30, 50, 10)
  
  Z %>%
  dplyr::mutate(coef = factor(coef, levels = c("Difference", "Overall mean"))) %>%
  ggplot(., aes(x = x, y = coef.mean)) +
  annotate("rect", ymin = -Inf, ymax = Inf, xmin = 35, xmax = 50, fill = "steelblue", alpha = 0.2) +
  geom_ribbon(aes(ymin = lower.mean, ymax = upper.mean), alpha = 0.2) +
  geom_path(linewidth = 1.2) +
  geom_hline(aes(yintercept = 0), color = "grey", linetype = "dashed") +
  geom_vline(aes(xintercept = 0), color = "grey", linetype = "dashed") +
  theme_bw() +
  facet_grid(coef ~ ., scales = "free_y") +
  theme(strip.text = element_text(face = "bold", size = rel(1.0))) +
  scale_x_continuous(name = "Days relative to sowing", breaks = breaks, labels = labels) +
  theme(axis.title.x = element_text(face = "bold", size = 11)) +
  ylab("Coefficient function") +
  theme(axis.title.y = element_text(face = "bold", size = 11))
  }

Function-on-scalar regression

Dew point

fANOVA(dat = data.prep(d2m))

using seWithMean for  s(yindex.vec) .

Mean air temperature

fANOVA(dat = data.prep(t2m_mean))

using seWithMean for  s(yindex.vec) .

Max air temperature

fANOVA(dat = data.prep(t2m_max))

using seWithMean for  s(yindex.vec) .

Min air temperature

fANOVA(dat = data.prep(t2m_min))

using seWithMean for  s(yindex.vec) .

Max - Min air temperature

fANOVA(dat = data.prep(tdiff))

using seWithMean for  s(yindex.vec) .

Temperature-Dewpoint depression

fANOVA(dat = data.prep(dpd))

using seWithMean for  s(yindex.vec) .

Soil temperature

fANOVA(dat = data.prep(st))

using seWithMean for  s(yindex.vec) .

Soil moisture

fANOVA(dat = data.prep(sm))

using seWithMean for  s(yindex.vec) .

soil temperature:soil moisture ratio

fANOVA(dat = data.prep(stsm))

using seWithMean for  s(yindex.vec) .

Surface pressure

fANOVA(dat = data.prep(sp))

using seWithMean for  s(yindex.vec) .

Relative humidity

fANOVA(dat = data.prep(rh))

using seWithMean for  s(yindex.vec) .

VPD

fANOVA(dat = data.prep(vpd))

using seWithMean for  s(yindex.vec) .

Growing degree days

fANOVA(dat = data.prep(gdd))

using seWithMean for  s(yindex.vec) .

# We will illustrate fdatest with the d2m variable.

# For fdatest, we need to create two separate matrices (for wm = 0 and wm = 1), which we can then put into a list for convenience.
wm0 <- 
  wm_data %>%
  dplyr::filter(wm == -1) %>%
  dplyr::select(subject, dap, d2m) %>%
  tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = d2m) %>%
  dplyr::select(-subject) %>%
  as.matrix()
# and ...
colnames(wm0) <- NULL

wm1 <- 
  wm_data %>%
  dplyr::filter(wm == 1) %>%
  dplyr::select(subject, dap, d2m) %>%
  tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = d2m) %>%
  dplyr::select(-subject) %>%
  as.matrix()
colnames(wm1) <- NULL

# Performing the ITP:
ITP.result <- fdatest::ITP2bspline(wm0, wm1, B = 100)

# The function generates a print line for each iteration. To suppress that, wrap within sink:
tic()
{ sink(type = "message"); ITP.result <- fdatest::ITP2bspline(wm0, wm1, B = 1000); sink() }
toc()  # 9.18 sec

# Plotting the results of the ITP:
# (there are two plots. The first is of the individual curves. The 2nd is of the adjusted p-values)
plot(ITP.result, main = NULL, xrange = c(1, 365), xlab = 'Day')

# Plotting the p-values heatmap
# (I'm not finding it all that telling)
ITPimage(ITP.result, abscissa.range = c(0, 1))

# Selecting the significant components at 5% level:
which(ITP.result$corrected.pval < 0.05)

# Which corresponds to the following days (relative to sowing, where sowing is day = 0):
seq(-30, 50)[31:32]

# NEXT:
# Take the above code, and wrap into functions to (i) prep the data for input into fdatest, (ii) calling the fdatest functions to estimate the adjusted p values

make_kable <- function(...) {
  # kable and kableExtra styling to avoid repetitively calling the styling over and over again
  # See: https://stackoverflow.com/questions/73718600/option-to-specify-default-kableextra-styling-in-rmarkdown
  # knitr::kable(...) %>%
  kable(..., format = "html", row.names = TRUE, align = 'l') %>%
    kable_styling(bootstrap_options = c("striped"), position = "left", font_size = 11, full_width = FALSE) 
}

do.fdatest <- function(x) {
  # Performs an Interval Testing Procedure for testing the difference between the two functional wm populations evaluated on a uniform grid
  #
  # Args:
  #  x = unquoted variable name, e.g. d2m
  #
  # Returns:
  #   a tibble of the variable with a list vector of the days (relative to sowing) where the two populations differ functionally
  
  # Data prep for input to fdatest
  .x <- enquo(x)
  # wm0 = N x J matrix of the functional data for wm absent
  wm0 <- 
    wm_data %>%
    dplyr::filter(wm == -1) %>%
    dplyr::select(subject, dap, !!.x) %>%
    tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = !!.x) %>%
    dplyr::select(-subject) %>%
    as.matrix()
    # and ...
    colnames(wm0) <- NULL
    
  # wm1 = N x J matrix of the functional data for wm present  
  wm1 <- 
    wm_data %>%
    dplyr::filter(wm == 1) %>%
    dplyr::select(subject, dap, !!.x) %>%
    tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = !!.x) %>%
    dplyr::select(-subject) %>%
    as.matrix()
    colnames(wm1) <- NULL
  
  dat <- list(wm0 = wm0, wm1 = wm1)
  
  # The function generates a print line for each iteration. To suppress that, wrap within sink:
  # (set a seed for reproducibility)
  { sink(nullfile()); set.seed(86754309); foo <- ITP2bspline(dat$wm0, dat$wm1, B = 10000); sink() }

  # Selecting the significant components at 5% level:
  # Which corresponds to the following days (relative to sowing, where sowing is day = 0):
  z <- seq(-30, 50)[which(foo$corrected.pval < 0.05)]
  
  res <- tibble(var = rlang::as_name(.x), days = list(z))
  
  return(res)
}

# Example of use:
# do.fdatest(x = d2m)

filter.iwt <- function(x) {
  # Filter the interval-wise testing results to see the days that were significant
  # Args:
  #  x = the series (quoted character string)
  # Returns:
  #  a vector of the days where the series was different between wm = 0 and wm = 1
  iwt |>
  dplyr::filter(var == x) |>
  purrr::pluck("days", 1)
}


zee <- function(series) {
  # Output the start and end days of significant windows within a time series
  # Args:
  #  series = quoted character name of the series
  # Returns:
  #  a table
  v <- 
    filter.iwt(x = series) %>% 
    split(., cumsum(c(1, diff(.) != 1)))
  
  # Create an empty data frame with three named columns:
  z <- data.frame(matrix(
    vector(), 0, 3, dimnames = list(c(), c("series", "start", "end"))), 
    stringsAsFactors = F)
  
  # Now loop over v to pick out the start and end of the continuous windows:
  for (i in 1:length(v)) {
    b <- purrr::pluck(v, i) %>%  dplyr::first()
    c <- purrr::pluck(v, i) %>%  dplyr::last()
    z[i, "series"] <- series
    z[i, "start"] <- b
    z[i, "end"] <- c
    } # end for loop
  return(z)
}

# Examples of use:
# zee("t2m_mean")
# zee("d2m")

#  Example of the workflow:
###---###
# tst.d2m <- do.fdatest(x = d2m)

# You could output this way:
# pluck(iwt.list, "tst.d2m") %>% 
#   make_kable()

# But this shows the windows in a cleaner format:  
# zee("d2m") %>% 
  # make_kable()
###---###

# Now we're ready to roll...

load(here::here("FunctionalDataAnalysis", "FunctiononScalar", "iwt.RData"))

Interval-wise tests

Perform the interval-wise tests

tst.d2m <- do.fdatest(x = d2m)

tst.t2m_mean <- do.fdatest(x = t2m_mean)

tst.t2m_max <- do.fdatest(x = t2m_max)

tst.t2m_min <- do.fdatest(x = t2m_min)

tst.tdiff <- do.fdatest(x = tdiff)

tst.dpd <- do.fdatest(x = dpd)

tst.st <- do.fdatest(x = st)

tst.sm <- do.fdatest(x = sm)

tst.stsm <- do.fdatest(x = stsm)

tst.sp <- do.fdatest(x = sp)

tst.rh <- do.fdatest(x = rh)

tst.vpd <- do.fdatest(x = vpd)

tst.gdd <- do.fdatest(x = gdd)

I’ve already done the interval-wise tests, so the chunks in this section will just output the days that the test found to be significant.